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BOUNDARY  AND  INTERFACE  CONDITIONS  FOR  HIGH  ORDER  FINITE 
DIFFERENCE  METHODS  APPLIED  TO  THE  EULER  AND  NAVIER-STOKES 

EQUATIONS 

JAN  NORDSTROM  *  AND  MARK  H.  CARPENTER  t 

Abstract.  Boundary  and  interface  conditions  for  high  order  finite  difference  methods  applied  to  the 
constant  coefficient  Euler  and  Navier-Stokes  equations  are  derived.  The  boundary  conditions  lead  to  strict 
and  strong  stability.  The  interface  conditions  are  stable  and  conservative  even  if  the  finite  difference  operators 
and  mesh  sizes  vary  from  domain  to  domain.  Numerical  experiments  show  that  the  new  conditions  also  lead 
to  good  results  for  the  corresponding  nonlinear  problems. 

Key  words,  high-order  finite-difference,  numerical  stability,  interface  conditions,  summation-by-parts 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  In  many  computational  problems,  low  order  finite  difference  methods  (second  order 
or  less)  axe  not  accurate  enough.  Examples  in  which  high-frequency  components  in  the  solution  must  be 
resolved  by  using  high  order  finite  difference  methods  (HOFDM)  include  aeroacoustics,  turbulence  and 
transition  simulations,  the  propagation  and  scattering  of  electromagnetic  waves,  and  simulation  of  reactive 
flows  at  high  speeds  [1],  [2],  [3],  [4],  [5],  [6].  The  efficiency  [7]  of  HOFDM  can  be  used  either  to  increase 
the  accuracy  for  a  fixed  number  of  mesh  points  or  to  reduce  the  computational  cost  for  a  given  accuracy  by 
reducing  the  number  of  mesh  points. 

The  main  reason  that  low  order  finite  difference  methods  are  used  in  practical  calculations  is  because  of 
the  difficulty  that  arises  for  HOFDM  near  the  boundaries  of  the  computational  domain.  On  a  Cartesian  mesh, 
it  is  quite  easy  to  derive  nonsymmetric  boundary  operators  that  have  high  formal  accuracy,  the  difficulty  is 
to  derive  highly  accurate  and  stable  operators.  In  [8]  and  [9]  HOFDM  are  constructed  based  on  the  work  in 
[10]  and  [11].  In  these  strictly  stable  schemes,  the  growth  rates  of  the  analytic  and  semidiscrete  solution  are 
identical.  Strict  stability  is  obtained  by  constructing  discrete  operators  that  satisfy  a  summation-by-parts 
(SBP)  rule  which  mimics  the  integration-by-parts  rule  in  the  continuous  case.  For  calculations  over  long 
times,  strict  stability  is  very  important  beacuse  it  prevents  error  growth  in  time  for  fixed  Ax. 

In  [12]  it  was  shown  that  many  G-K-S  stable  [13]  (convergence  to  true  true  solution  as  Ax  0)  scalar 
schemes  were  not  strictly  stable.  Moreover,  many  scalar  schemes  that  were  both  G-K-S  stable  and  strictly 
stable  exhibit  time  growth  when  they  are  applied  to  systems  of  equations.  The  underlying  reason  for  the 
error  growth  in  time  caused  by  the  way  the  mathematical  boundary  conditions  were  imposed.  An  orthogonal 
projection  operator  is  used  to  impose  the  mathematical  boundary  conditions  in  [14]  and  [15].  In  the  so  called 
SAT  (simultaneous  approximation  term)  procedure  [16],  a  linear  combination  of  the  boundary  conditions  in 
the  form  of  a  forcing  function  and  the  differential  equations  is  solved  near  the  boundary.  Both  these  methods 
impose  the  correct  boundary  conditions  and  preserve  strict  stability. 
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Administration  under  NASA  Contract  No.  NASl-97046  while  this  author  was  visiting  the  Institute  for  Computer  Applications 
in  Science  and  Engineering  (ICASE),  Mail  Stop  403,  NASA  Langley  Research  Center,  Hampton,  VA  23681-2199. 
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Another  important  concept  is  strong  stability.  An  approximation  is  strongly  stable  if  the  solution, 
including  the  values  at  the  boundary  points,  can  be  estimated  in  terms  of  all  data  in  the  problem  [17],  The 
stability  estimate  in  the  strongly  stable  case  leads  directly  to  the  error  estimate  if  no  extra  or  numerical 
boundary  conditions  are  necessary.  Stability  analysis  using  the  Laplace  transform  technique  leads  to  strong 
stability  if  the  Kreiss  condition  is  satisfied;  see  [18]  and  [9]  paper  IV.  Note  that  strict  stability  leads  to  strong 
stability,  but  strong  stability  does  not  imply  strict  stability. 

Most  investigations  regarding  HOFDM  are  done  on  linear  hyperbolic  model  equations  with  constant 
coefficients  on  a  uniform  mesh.  However,  nonlinear  Navier-Stokes  calculations  on  nonuniform  meshes  have 
been  performed  [19].  One  of  the  conclusions  in  [19]  was  that  the  treatment  of  the  metric  derivatives  is  a 
crucial  point  for  nonsmooth  meshes.  This  problem  is  analyzed  in  [20]  where  so  called  mimetic  difference 
operators  (discrete  operators  with  the  same  symmetry  properties  as  the  continuous  operators)  are  derived. 
In  [15],  strict  stability  for  parabolic  and  hyperbolic  systems  in  curvilinear  coordinates  on  a  single  domain 
were  investigated. 

Generating  a  smooth  grid  around  a  complex  configuration  can  be  very  difficult,  if  not  impossible,  and  is 
often  the  most  time-consuming  aspect  of  the  solution  procedure.  This  fact  has  limited  the  use  of  HOFDM  in 
practical  calculations  to  the  small  class  of  simple  geometries  which  can  be  smoothly  mapped  onto  the  unit 
cube.  In  this  paper  we  consider  a  structured  multiblock  approach  in  which  each  subdomain  is  discretized  by 
using  a  discrete  operator  with  the  SBP  property.  The  subdomains  are  patched  together  to  a  global  domain 
by  using  suitable  interface  conditions.  This  technique  was  used  in  [21],  [22]  and  [23]  for  Chebyshev  spectral 
methods. 

In  [24],  stable  and  conservative  interface  conditions  for  HOFDM  applied  to  the  scalar  advection-diffusion 
equation  on  multiple  domains  were  derived.  In  each  subdomain  the  step  size  was  constant  but  significantly 
different  from  that  in  the  adjacent  subdomains.  Also,  the  finite  difference  operators  could  vary  from  sub- 
domain  to  subdomain.  In  this  paper  we  will  generalize  the  results  in  [24]  and  extend  the  analysis  to  the 
one-dimensional  constant  coefficient  Euler  and  Navier-Stokes  equations. 

The  rest  of  this  paper  will  proceed  as  follows.  In  section  2,  some  basic  definitions  are  given.  In  section  3, 
the  Navier-Stokes  equations  on  conservative,  primitive,  and  characteristic  variable  form  are  given.  In  section 
4,  the  continuous  problem  is  analyzed,  while  the  discrete  problem  is  investigated  in  section  5.  Numerical 
experiments  are  performed  in  Section  6  and  we  summarize  and  draw  conclusions  in  section  7. 

2.  Definitions.  Consider  the  linear  initial  boundary  value  problem 


Wt 

=  Pw  -{-5F{x,t) 

,x  € 

,t>0, 

w 

=  Sf{x) 

,x  e 

,t  =  o, 

Lcu) 

11 

,x  e  r 

O 
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where  P  is  the  differential  operator  and  Lc  is  the  boundary  operator.  The  initial  function  6f,  the  forcing 
function  SF,  and  the  boundary  data  5g  are  the  data  of  the  problem;  w  denotes  the  difference  between  a 
solution  with  data  f^F^g  and  one  with  data  /  -h  (5/,E  +  J/ ^g-^Sg.  There  are  many  concepts  of  well  posedness, 
see  [17].  Here  we  consider  the  following  definition. 

Definition  1.  The  problem  (2.1)  is  strongly  well  posed  if  the  solution  w  is  unique,  exists,  and  satisfies 

(2.2)  \\w\\l  +  f  Midi  <  +  \\Sg\\l)dt}, 

Jo  Jo 

where  Kc  and  rjc  may  not  depend  on  5F,  5f,5g.  ||  •  \\q  and  ||  *  ||r  suitable  continuous  norms. 
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The  semidiscrete  version  of  (2.1)  is 

{wj)t  —  Qwj  -h  SFj{t)  ^Xj  gQ,  >  0, 

(2.3)  Wj  =  5fj-  ,XjGFt  ,t  =  0, , 

Ld'Wj  —  5g{t)  ^Xj  gT  ,t>0, 

where  Q  is  the  difference  operator  approximating  the  differential  operator  P,  5Fj  is  the  forcing  function,  5fj 
the  initial  function,  Ld  the  discrete  boundary  operator  where  numerical  boundary  conditions  are  included, 
and  Sg  the  boundary  data.  It  is  assumed  that  (2.3)  is  a  consistent  approximation  of  (2.1). 

Closely  related  to  the  concept  of  well  posedness  is  the  concept  of  stability. 

Definition  2.  The  problem  (2.3)  is  strongly  stable,  if  for  a  sufficiently  fine  mesh,  the  solution  wj 
satisfies 

(2.4)  ||«;||^  +  f  Midi  <  +  f\\\5F\\l  + 

Jo  Jo 

where  Ka  and  rjd  may  not  depend  on  5Fj,  Sfj,Sg.  ||  •  \\q  and  1|  •  ||r  are  suitable  discrete  norms. 

Definition  3.  The  approximation  (2.3)  of  (2.1)  is  strictly  stable  if  the  analytical  and  discrete  growth 
rates  (see  (2.2)  and  (2.4))  satisfy 

(2.5)  7]d<Vc  +  0{Ax), 
where  Ax  is  the  mesh  size. 

For  later  reference  we  also  define  some  useful  matrix  operations;  see  [25]. 

Definition  4.  Let  A  be  a  p  x  q  matrix,  B  be  an  mxn  matrix,  and  Ii  the  I  x  I  identity  matrix,  then 


A®B  = 

^  Uo,oP 

ao,g  iB  ^ 

,  Ii<S>  B  — 

/  B  0  ■ 

0  B  ■ 

■  ^  \ 
■  0 

\  ap  i,oP 

‘  *  ap  iB  j 

•  o 

•  o 

■  B  ) 

The pxq  block  matrix  A<^B  and  the  I  xl  block  diagonal  matrix  // are  called  Kronecker products.  There 
are  a  number  of  rules  for  Kronecker  products  (see  [25]).  In  this  paper  we  will  make  use  of, 

(2.6)  {A®B){C®D)  =  {AC)^iBD),  {A®  Bf  =  A^  ®  B'^ . 

The  following  lemma  will  be  used  frequently  below;  it  is  a  direct  consequence  of  the  first  rule  in  (2.6). 

Lemma  1.  Let  A  be  anmxm  matrix,  B  be  an  nxn  matrix,  A  =  In^A  and  B  =  B^Im,  then  AB  =  BA. 
Proof  :  The  first  condition  in  (2.6)  leads  to  AB  —  {In  0  A)(P  0  Im)  =  P  0  A  =  (P  0  Im){In  0  A)  =  BA. 
□ 

3.  The  Euler  and  Navier-Stokes  equations.  The  one-dimensional  constant  coefficient  Navier- 
Stokes  equations  in  primitive  (VF),  characteristic  (C),  and  conservative  (Q)  variable  form  are 

(3.1)  Wt-\-AW^  =  eBW^^,  Ct  +  AC^  =  eXC^^,  Qt  +  F^  =  €F^ 

respectively.  With  e  =  0,  equation  (3.1)  becomes  the  one-dimensional  constant  coefficient  Euler  equations. 
The  overbar  is  used  to  denote  variables  with  a  constant  state.  The  relation  between  W,  C,  Q  where  W  = 
(p,  u,  TY  is 

(3.2)  C  =  RSW,  Q  =  TW 
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where 
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/  1  0  0 
u  p  0 

\  cV(7(7  -  1))  +  «V2  pu  p/ilil  -  1)M^) 


Note  that  RR^  —  h- 

The  transformation  (3.2)  implies  that  the  matrices  and  fluxes  in  (3.1)  are 


(3.3) 


A  = 


fu  p  0  ^ 

^hP  U  1/7-^to 

\  0  (7  -  u  / 


(3.4) 


B 


/  0  0  ^  \ 
0  (A  +  2p,)/ p  0 
0  0  jRf{Prp)  j 


(3.5) 


a  —  c  0  0 


A  =  {RS)A{RS)  ^  =  I  0  u  0 


0 


u  +  c 


(3.6) 


X  =  {RS)B{RS)  1  =  ^ 


f  9  +  ^  6  -  ^  \ 

—acj) 

^  6  —  ^  ^  +  0  / 


(3.7) 


=  TAW  =  TAT  =  TBW^  =  TBT 


The  dependent  variables  and  parameters  p,u,T,p,c,  Moo,Ai,A,«;,Pr,7  and  e  are  respectively  the  density, 
x^y^z  components  of  the  velocity,  the  temperature,  the  pressure,  the  speed  of  sound,  the  free-stream  Mach 
number,  the  shear  and  second  viscosity,  the  coefficient  of  heat  conduction,  the  Prandtl  number,  the  ratio  of 
specific  heats,  and  the  inverse  Reynolds  number.  The  notations  0  =  (A  +  2p)fp,^  =  (7  “  = 

—  1)  has  also  been  introduced. 

4.  The  continuous  problem.  In  this  paper  we  will  consider  interface  conditions  between  subdomains. 
However,  interface  conditions  are  closely  related  to  boundary  conditions;  therefore,  we  start  with  the  single 
domain  problem. 
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4.1.  The  continuous  single  domain  problem.  To  make  the  presentation  self-contained,  some  results 
in  [27]  are  included  in  this  section.  Consider  the  Navier-Stokes  equations  on  characteristic  form, 

Ci  +  KCx  =  "I" 

,i  =  0  ,-l<x<l, 

,t>0  ,x  =  -l, 

,i>0  ,a;  =  -hl, 

where  C  —  {pcu  -  p,  a{pc^  -  p),  pcu  -h  p)^,0  <  e  «  1  and  L  i,  T+i  are  the  boundary  operators.  For  fZ  >  0, 
there  is  inflow  at  x  —  —1  and  outflow  at  a:  =  1. 


(4.1) 


C-  =  fix) 

L  iC  =  9  lit) 
L+iC  =  g+iit) 


4.1.1.  Well  posedness.  Let 

iu,  V)  =  U^Vdx,  iU,  U)  =  \\uf,  Ilt/||2  =  \u\l^  1  +  \u\U+i 

denote  the  L2  scalar  product,  the  L2  norm,  and  the  boundary  norm  respectively.  The  energy  method  applied 
to  (4.1)  leads  to 

lldl?  =  [C^AC  -  2eC^XCa,]l=^\  -  2e(C„  XC,)  +  2(C,  F). 


The  boundary  conditions  (see  [27]  and  [22]) 

(4.2) 


L  _  eXCx  =  5  1, 


L+iC  =  {^^  J^^^C-eXCA  ={9+i}i,  i  =  l,2, 


(4.3)  L+iC  =  |-  2 

where  |A|  =  diagUXil,  [Az],  [As])  leads  to 

IlCjjf  =  -2e(a,Xa)  +  2(C,F) 

(4.4)  —[C^AiC  —  2C^g  i]a:=  \  AqC  ■¥  2C'^ g+i\x=+ii 

where  g+i  =  (fifi, 92,51  -  (2/Q!)92)^  and 

jAij  0  (|Ai|-Ai)/2 


(4.5) 


A/  =  |A|,  Ao  = 


0  [Azl  0 

(|Ail-Ai)/2  0  jAsI 


Integration  of  (4.4)  leads  to 

\\Of  +  e’<'^{2e j\c„XC,)e  ’'‘*  +  5j(  IIC#e  ’“*} 

(4^6)  <  e-«’{ll/r  +  j  llsllfe  "‘df  +  1  jf''  ||f  fe  -*}, 

where  0  <  77  <  1,  <5  =  min|dj|,  D  =  |A|H,  and 

rr-^-  m  1  rri  rr  _  MzM  rr  =  MziM 
H  -dtagiHi,l,H3),  Hi  3  |;>,3|  +  |Ai|- 

Note  that  (4.2), (4. 3)  reduce  to  the  characteristic  boundary  conditions  for  the  Euler  equations  as  e  — ^  0. 

Uniqueness  follows  directly  from  the  estimate  (4.6).  Existence  can  be  shown  by  using  the  Laplace- 
transform  technique  or  via  difference  approximations;  see  [26]  and  [28].  Since  (4.6)  is  of  the  form  (2.2),  we 
can  conclude  that  the  following  theorem  holds. 

Theorem  1.  The  problem  (4-V  boundary  conditions  (4-2), (4’3)  i'S  strongly  well  posed. 
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4.2.  The  continuous  multiple  domain  problem.  In  this  section  we  split  the  domain  [-1,1]  into 
[—1,0]  and  [0, 1]  and  focus  on  the  interface  problem  at  a;  —  0.  The  two  coupled  problems  are 


(4.7) 


Ut  +  XU,  —  tXUxx  +  F(Xi  t) 

,t>0  ,-l<x<0, 

u  =  f(x) 

,t  =  0  ,-l<x<0, 

II 

>  0  ,x  =  -1, 

Lo(U-V)  =  0 

,t>0  ,x  =  0, 

Vt  +  AV,  =  eXV„  +  F(x,t) 

,<>0  ,0<x<+l, 

V  =  f(x) 

,t  =  0  ,0<x<+l, 

Lo(V-U)  =  0 

cs- 

IV 

0 

II 

0 

L+iV  =  g+i(t) 

>  0  ,x  =  +1, 

(4.8) 


respectively.  The  characteristic  variables  in  the  left  [-1, 0]  and  right  [0,  +1]  domain  are  U  and  V  respectively. 
The  coupling  between  (4.7)  and  (4.8)  is  given  by  the  operator  Lq. 

By  subtracting  (4.1)  from  (4.7-4.8),  by  transforming  the  problem  on  [0, +1]  onto  [—1,0]  via  the  trans¬ 
formation  X  — >  — and  finally  by  replacing  ^  with  x,  we  obtain 


(4.9) 

where 

and 

(4.10) 


+  A'i/jx  =  eXijjx 

Ip  =  0 

L  lU  =  0 
L+iV  =  0 
Lo{U-V)  =  0 


,t>0  ,-l<x<0, 
,t  =  0  ,  -1  <  a:  <  0, 
,t>0  ,x  =  -1, 

>  0  ,x  =  -1, 
,t>0  ,x  =  0, 


iP 


U 

V 


u-c 

v-c 


A  = 


+A  0 
0  -A 


X 


X  0 
0  X 


L  ^  \  U  -  eXU^,  L+iV 


I  = 

)  i 


2  -  - — TJ. '  j  2 

4.2.1.  Well  posedness.  The  energy  method  applied  to  (4.9)  leads  to 

=  [V'^AV'  -  2ei/)^XV>x]x=o  ^  -  2e(V’x,  Xip^). 

The  analysis  of  the  single  domain  problem  implies  that  the  boundary  terms  at  x  =  —  1  are  negative  semidef- 
inite  with  the  boundary  operators  (4.10).  At  the  interface  x  =  0,  we  have 


'ilFK'ij)  -  2e'ilF X\l)x]x^o  = 


(  u-v  \ 


(4.11) 


U  +  V 

(U-V), 

\  (U  +  V)a.  j 


( 


V 


0 

A 

-tX 

0 


A  -tX 
0  0 

0  0 

-eX  0 


0  \ 


(  U-V  \ 


U  +  V 
(U  -  V), 

V  (U  +  V),  ) 


-eX 
0 

0  / 

Well  posedness  for  the  Euler  equations  (e  =  0)  requires  {7  —  ^  =  0  since  A  is  nonsingular.  With  that 
choice  we  get 

[ip'^K^  -  2€iFX^,],=Q  =  -2eU'^X(U  +  V),  =  -2e((RS)'^U)^ B(Wi  +  W2),, 
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where  {RS)  ~  Wl,  {RS)  —  Wr  denotes  the  primitive  variables  in  the  left  and  right  domain  respec¬ 
tively.  The  structure  of  B  (see  (3.4))  and  a  transformation  to  the  original  coordinate  system  lead  to  the 
following  theorem. 

Theorem  2.  If  Theorem  1  holds  and  the  interface  conditions 


(4.12) 


U-V 
{U  -  v% 


1  0  1  \ 

1  a  -1  y 


are  used,  then  (4-7)  and  (4-8)  are  strongly  well  posed. 

Remark.  The  problems  (4.7)  and  (4.8)  axe  strongly  well  posed  in  the  sense  that  the  solutions  can  be 
estimated  in  ternms  of  the  data  in  the  corresponding  one  domain  problem  (4.1). 

Remark.  The  condition  (4.12)  in  primitive  variable  formulation  is 

{  Uo, 

\tD,  J\  {Wl  -  IVr),  )  \0  f>  1  J 


5.  The  discrete  problem.  Let  UfDU  be  the  numerical  approximations  of  the  scalar  quantities  u  and 
Ux  respectively.  The  approximation  VU  of  the  first  derivative 

VU  =  P  ^QU,  Pu^-Qu  =  PTei,  |Tei|  =  e)(Ax™,Ax") 


satisfies  the  SBP  rule 

(5.1)  {U,  VV)p  =  UnVn  -  UoVo  -  {VU,  V)p 
where 

(5.2)  iU,V)p  =  U^PV,  P  =  P^,  Q  +  Q'^  =  D,  D  =  diag[-l,0,  ..,0,1] 

and  0  <  pminl^xl  <  R  ^  pmaxI^xT  Operators  of  the  SBP  type  arise  naturally  with  centered  difference 
approximations;  for  examples  see  [11], [29],  [12], [30]. 

The  second  derivative  can  be  obtained  by  applying  the  first  derivative  operator  twice.  Such  an  ap¬ 
proximation  satisfies  the  SBP  rule  (5.1)  exactly.  However,  there  are  drawbacks  with  such  a  procedure.  A 
second  derivative  formed  in  that  way  is  unnecessarily  wide  and  inaccurate  and  can  lead  to  odd-even  mode 
decoupling.  A  second  derivative  operator  with  the  following  properties, 

(5.3)  V‘^U  =  P  ^RU,  Pu^^-Ru  =  PTei,  Tea  -  0(Ax'",  Ax"), 

(5.4)  R  =  (-5^M  +  D)S, 

was  suggested  in  [24].  The  matrix  D  is  given  in  (5.2);  M  is  positive  definite,  i.e.,  U^MU  >  0  and  0  < 
rUminAx/  <  M  <  mmax^xl. 

5  is  a  diagonal  matrix  with  a  discrete  representation  of  the  first  derivative  on  the  first  and  last  rows, 

~  ^a;(^0jt)  Te3,  —  V,x{^n^t^  “1“ 
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where  ITesI  =  C’(Aa;’')  and 


SOO  SOl  S02  S03 

0  1  0 

0  1  0 

0  10 
0  10 

*  ’  ’  ^nn  3  ^nn  2  ^nn  1  ^nn 

The  second  derivative  defined  in  (5.3)  and  (5.4)  satisfies  a  modified  SBP  rule  We  have 

(U,  P V)p  =  Un{VV}n  -  Uo{VV}o  -  {SUfM{SV). 

The  notation  |Tei|,  |Te2|  =  0{Ax'^,Ax'^)  and  iTesl  “  0{Ax'^)  means  that  the  approximation  of  the 
differential  operator  is  accurate  to  order  m  in  the  interior  of  the  domain,  to  order  n  at  the  boundary  and 
that  the  approximation  of  the  boundary  conditions  is  accurate  to  order  r.  The  relation  between  the  different 
orders  of  accuracy,  i.e.,  m,  n,  r  is  discussed  in  section  5.1.2  below. 

Examples  of  first  and  second  derrivative  approximations  are  given  in  (A.1)-(A.5)  in  appendix  A.  The 
approximations  are  second  order  accurate  in  the  interiour  of  the  domain  and  first  order  accurate  at  the 
boundary.  This  means  that  for  (A.1)-(A.5)  we  have  m  =  2  and  n  =  1. 

So  far  we  have  considered  difference  approximations  of  scalar  quantities.  The  corresponding  approxima¬ 
tions  for  vector  quantities  are  defined  by  using  Kronecker  products  (see  definition  4).  The  spatial  operators 
X),  and  the  matrices  that  define  them  are  of  the  form  B  0  is  in  this  paper.  As  an  example,  P  means 
{P  ^  0  /3)(Q  0  /a)  =  P  0  Is.  In  the  sequel,  that  notation  is  implied. 

Let  H  —  >  0;  for  later  reference  we  introduce  the  notations 

(5.5)  {U,V)h  =  U^HV,  {U,U)h^\\U\\%,  \\U\\l^  = 

5.1.  The  discrete  single  domain  problem.  We  introduce  a  uniform  mesh  Xi  =  —1  H-iAx,xo  — 
—  lyXn  =  +1-  The  finite  difference  approximation  of  (4.1)  with  the  SAT  technique  [16]  for  boundary  condi¬ 
tions  is 

Ct-\-AVC  =  eXV^C^F 

(5.6)  +  P  ^{a  i{L^,C-g  i)e  i  +  a+i{L^,C  -  g+i)e+i}, 

C(0)  =  f, 

where 

(5.7)  V  =  P  ^Q®h, 

(5.8)  =  or  R=  {-S'^M  +  D)S  ®  h, 

(5.9)  A  =  /„ig)A,  X  =  In<E)X,  e  i  =  (1,  ...0)^  ® /a,  e+i  =  (0,  ...1)^  (8»  Za. 

The  unknown  diagonal  matrices  a  i  and  cr+i  will  be  determined  below. 


5.1.1.  Stability.  The  energy  method  leads  to 

^IICII^  =  -C'^ikQ  +  Q'^l)C  +  eC'^iXR  +  B?'X)C  +  2{C,  F)p 
at 

(5.10)  -\-2Cq(T  i\L^]C  —  g  i]  +  2C'^cr+i[i^iC  —  54.1]. 

The  definition  of  the  first  derivative  operator  P  and  Lemma  1  leads  to 

(5.11)  -C'^iAQ  +  Q^A)C  =  C'^AC'o  -  C^ACn- 

The  definition  of  the  second  derivative  operators  R  =  {-S^M  +  D)S  and  R  =  QP  yields 

C'^{XR  +  R^X)C  =  -2CqXVCo  +  2CnXVCn 

(5.12)  -  {SCfiXM  +  {XMf){SC) 

C'^iXR  +  R'^X)C  =  -2C0XVC0  +  2CnXVCn 

(5.13)  -2{P  ^QCYPXP  ^QC, 

respectively. 

By  introducing  (5. 11), (5. 12)  and  (5.13)  into  (5.10)  we  get 

^IIC^IIp  +  2e{VC,  XVC)h  =  [C^AC  -  2tC'^XVC^=l  +  2{C,  F)p 
at 

+  2C^a  i[LV-5  i] 

(5.14)  +  2Cj;cr+i[L^lC  —  5+1], 
where  the  scalar  products  and  norms  are  defined  in  (5.5)  and 

R  =  QP  ^Q=^H  =  P, 

R={-S^M  +  D)S^H  =  {SiP  ^Q)  Yi^\^)iS{P  ^Q)  ^)- 

The  boundary  operators  are  the  discrete  versions  of  (4.2)-(4.3),  with  one  important  modification. 

In  [27]  it  is  shown  that  the  two  outflow  conditions  in  (4.3)  determine  the  value  of  the  last  row  of  XCx  in 
terms  of  the  in-going  characteristic  variable  and  boundary  data;  i.e.,  (4.3)  implies  that 

(5.15)  {-em}3  =  -^i^^Ci+5r-(2/a)52,  x  =  +l. 

To  explicitly  incorporate  (5.15)  into  (5.6)  we  use 

(5.16)  -  eXVC^  =  g  i, 

(5.17)  L^,C  =  I  -  eXVC^  =  ff+i, 

where  (5'+i)3  is  equal  to  the  right-handside  of  (5.15).  The  boundary  conditions  (5. 16), (5. 17)  inserted  in 
(5.14)  yields 

^lldl^  =  -2e{VC,XVC)H  +  2{C,F)p 

-j-  cr  i(A  -}-  |A|)](7}i=o  +  {C^\—2eX  —  2e(j  iX\DC}iz=^Q 
+  {C'^[-A  +  a+i(A  -  |A|)]C}^^„  +  {C^[+2eX  ^  2ea^iX]VC}i=n 

(5.18)  +  {cr^iC'^C^(Ai  —  |Ai|)}i=n  -h  2CQg  i  -  2C^g^i. 
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The  choice, 


(5.19) 


a  I  —  -Is,  CT-f-i  =  Is, 


leads  to 


(5.20) 


lie'll?  =  -  2e{VC,XVC)H  +  2{C,F)p 

-  [C^AiC  -  2C^g  i]i=o  -  [C^fioC  +  2C^5+i]i=„, 


i.e.,  a  growth  rate  which  is  exactly  the  same  as  in  the  continuous  case  (compare  (5.20)  with  (4.4)).  The 
definitions  of  A/,  Aq  are  given  in  (4.5).  Integration  of  (5.20)  leads  to 

IIC^IIp  +  e^‘^^{2e  j\vC,XVC)He  \\C\\l^e 

(5.21)  <e’'"^{||/||?>  +  -^  r  Mhe  +  ^  T  \\F\\%e  ^-*dt}. 

Jo  rjD  Jo 

The  estimate  (5.21)  is  similar  to  (2.4)  and  hence  (5.6)  is  a  strongly  stable  approximation.  The  problem 
(5.6)  is  also  strictly  stable  (we  can  choose  tjd  ~  V  and  So  =  S,  see  (4.6), (2. 5)).  We  can  summarize  the  result 
in  the  following  way. 

Theorem  3.  The  approximation  (5,6)  of  the  problem  (4-i)  is  both  strictly  and  strongly  stable  if  (5.19) 
holds. 


5.1.2.  Accuracy.  The  problem  describing  the  deviation  Ej  ~  C{xj.,t)  —  Cj{t)  between  the  exact 
continuous  solution  and  the  discrete  approximation  given  by  (5.6)  is 


Et+kVE  =  eXV^E-\-T 

(5.22)  +  P  i{L^iE)e  i  +  a+i{L^,E)e+r}, 

E{0)  =  0. 

T  =  is  the  truncation  error.  and  comes  from  the  approximation  of  the  differential 

operator  and  the  approximation  of  the  boundary  conditions  respectively.  The  truncation  errors  have  the 
general  structure 


(5.23) 


rpDO 


0{Xx^)  \ 
e»(Ax'") 

e>(Ax'") 

0(Ax”) 


(  0(Ax(’"  1))  \ 
0 


nBC 


\  0{Ax^ 


(r  1) 


In  [31]  and  [32]  it  is  shown  that  difference  approximations  to  mixed  hyperbolic-parabolic  equations 
retain  the  accuracy  of  the  interior  scheme  {0{Ax'^))  if  a  finite  number  of  points  (independent  of  the  total 
number)  are  closed  with  boundary  stencils  {0{Ax^))  that  are  one  order  less  accurate.  A  requirement  for 
that  conclusion  is  that  an  energy  estimate  holds,  which  in  turn  means  that  the  mathematical  boundary 
conditions  must  be  approximated  to  the  order  of  the  internal  scheme.  The  discussion  above  implies  that 
that  n  —  m  —  1  and  r  =  m  is  necessary. 

We  will  now  apply  the  theory  in  [31]  and  [32]  to  the  type  of  difference  approximations  considered  in  this 
paper,  i.e.,  where  difference  operators  of  the  SBP  type  are  used  together  with  a  penalty  formulation  for  the 
boundary  conditions. 
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First,  we  split  E  and  the  T  into  two  parts,  i.e.,  E  =  +  E^  and  T  =  T^  +T‘^  where 


/ 


(5.24) 


= 


0 

9i 


\ 


9n  1 
0 


/  50 
0 


\ 


=  0{Ax^),  T2  = 


0 

\  9n  J 


=  C»(Aa;(’" 


Next,  we  use  the  energy  method  to  estimate  E'^.  The  energy  method  applied  to  (5.22)  with  E,  T  replaced 
by  E^,T^,  and  the  conditions  (5.19)  leads  directly  to 

IIE^IIp  <  0{Ax"^). 

Finally  we  use  the  Laplace- transform  technique  to  take  care  of  the  boundary  error  and  estimate  E^. 
So  far,  the  treatment  has  been  general.  However,  in  order  to  keep  the  algebraic  complexity  at  a  reasonable 
level,  we  now  need  to  simplify  and  be  specific.  We  will  consider  the  inviscid  (e  =  0)  Euler  equations  at  an 
inflow  boundary  approximated  with  the  second  order  scheme  given  by  (A.l)  and  (A.2)  in  appendix  A.  The 
half-plane  problem  obtained  by  Laplace-transforming  (5.22)  with  E,T  replaced  by  E^,T^  becomes 

sEl  +  A{El  -El)  =  a  1  (A -f  |A|)E2  +  Axg^, 

(5.25)  ~sEl  +  A(E?+i  -  E?  i)/2  =  0,  j  >  1, 

E]  ^  0,  j 

where  s  =  sAx.  The  second  and  third  equations  in  (5.25)  lead  to 


00 


(5.26) 


E] 


1^1  + cr2 


1^2  +  *^3 


/  0 
0 
VI 


^3’ 


(5.27)  i^j=  0  ,Aj  =  0  , 

-(s/Ai)  -  ^/TT(f/^  ,  Aj  <  0 

where  the  branch  of  the  square  root  is  the  one  with  positive  real  part  for  Re{s)  >  0.  The  case  when  Xj  —  0 
presents  no  problem;  it  only  reduces  the  number  of  equations  in  (5.6).  In  the  sequel,  we  assume  Xj  ^  0. 

The  first  equation  in  (5.25)  leads  to 

(5.28)  E{s)a  =  Axgo,  E{s)  =  diag{s  XjKj  -  {Xj  +  cr^  i{Xj  +  |Aj|)),  j  =  1,2,3. 

A  nonsingular  £7(5)),  i.e., 

(5.29)  det{E{s))  ^  0,  Re{s)  >  0, 
and  (5.26), (5. 27)  lead  to 

\Ej\  <  const. {Axgoly  i?e(s)  >  0,  j  >  0; 

i.e.,  the  Kreiss  condition  is  satisfied.  Parsevars  relation  and  the  fact  that  Ej{t)  cannot  depend  on  goiT)  for 
t  <T  leads  to 

f  <  const,  f  \Axgo\‘^dt^  j  >0 

Jo  Jo 
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Xfi  —  Xq 


axl 

^XR 

1  1 

) 

t  1.. 

1  1  .  I- 

- 1 - 1 - 1  ) — 

Xn  5 

- 1  1 — 

Xn  3 

Xn  ll 

1 

Xi 

X2  X^ 

X  =  0 


Fig.  5.1.  The  mesh  close  to  the  interface  atx  =  0. 


and  finally,  since  “  0(Ax), 


WE^Wp^OiAx^^). 

It  still  must  be  shown  that  (5.29)  holds.  The  inviscid  condition  for  strict  stability  Xj  ^{Xj  +  \Xj\)  <  0 
(see  (5.18)  and  (5.27))  which  implies  s+XjKj  =  lAjl^/l  +  {s/XjY  >  0,  leads  directly  to  (5.29).  The  procedure 
to  estimate  the  boundary  error  at  an  outflow  boundary  is  exactly  the  same  as  in  the  inflow  case.  We  can 
summarize  the  result  in  the  following  Theorem. 

Theorem  4.  The  approximation  (5.6)  of  the  problem  (4^1)  with  e  =  0  is  second  order  accurate  if 
Theorem  3  holds  and  the  first  derivative  operator  V  =  P  is  given  by  (A.l)  and  (A. 2)  in  appendix  A. 

Remark.  The  procedure  that  was  exemplified  above  to  prove  accuracy  in  the  second  order  accurate  case 
is  general.  The  last  step  where  one  uses  the  Laplace-transform  technique  to  estimate  the  boundary  error 

is  not  necessary  ^)  if  the  boundary  stencils  have  the  same  order  of  accuracy  as  the  internal  stencil,  i.e., 
n  =  m  and  ii)  if  the  approximation  of  the  mathematical  boundary  conditions  is  one  order  more  accurate, 
i.e.,  r  =  m  +  1. 


5.2.  The  discrete  multiple  domain  problem.  A  finite  difference  approximation  of  the  coupled 
problems  (4.7)  and  (4.8)  is 


Ut  +  KVlU 

= 

eXVlU  +  F  +  BTo 

+ 

Pl  -  Vo)  +  <yl{{VLU)n  -  {VRV)o))eL 

U{0) 

= 

f 

Vt  +  AVnV 

= 

eXVjiV  +  F  +  BTm 

+ 

PrH^W  -  Un)  +  <^U(PrV)o  -  {VLU)n))eR 

y(o) 

/■ 

The  characteristic  variables  in  the  left  (subscript  L)  [xq  =  -l,Xn  =  0]  and  right  (subscript  R)  [xq  =  0,Xm  = 
-hi]  domains  are  U  and  V  respectively,  see  Figure  5.1.  BTo.BTm  denote  the  boundary  terms  at  x  =  ±1 
respectively.  Definitions  of  X,e  i,e^i  are  given  in  (5. 7), (5. 8), (5. 9),  and  ep  =  (0, ,,  1)^  0  = 

The  values  of  cr  i  and  cr+i  that  lead  to  strict  and  strong  stability  for  the  discrete  single  domain 
problem  are  given  in  (5.19).  We  must  still  determine  Note  that  the  difference  operators 

Vl.V^.Vr.V^  can  be  different  in  the  left  and  right  domains  and  that  Axl  ^  Axr. 


5.2.1.  Conservation.  To  calculate  the  strength  and  speed  of  a  shock  with  finite  mesh  size,  one  needs 
a  conservative  scheme.  Let  us  start  by  considering  a  continuous  problem  in  conservation  form,  UtA'fx  = 
0,|x|<l,i>0.  Integration  over  the  domain  leads  to 

d 

—  y  ^  dx  +  Ui-  f  1=0, 
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i.e.,  the  total  change  of  u  in  the  domain  is  only  due  to  the  flux  through  the  boundaries.  Note  that  integration 
of  fx  over  the  the  domain  reverses  the  differentiation  process  and  leaves  information  only  at  the  boundaries. 
Let  F,  VF  denote  the  numerical  approximations  of  /,  The  discrete  SBP  derivative  satisfies 

(5.31)  fx-Vf  =  Teu  Vf  =  P^Qf,  rei  =  0(Ax"). 

Multiplying  (5.31)  with  the  operator  I'^P  where  =  [1, 1, ,  1]  (8)  /p  (/  has  p  components)  and  observing 
that  f^i  -  f  1  =  f^dx  leads  to 

i^Pfx  =  J  fx  dx  +  0{Ax'^). 

The  operator  l^P  is  the  discrete  integration  operator.  This  operator  reverses  the  process  of  differentiation, 
leaves  information  only  at  the  boundaries,  and  converges  to  the  continuous  integration  operator  as  Ax  ^  0. 
We  can  now  prove  the  following  theorem. 

Theorem  5.  The  approximation  (5.30)  of  the  problem  (4^1)  is  conservative  if 

(5.32)  (Tj^  —  —  A  =  0,  +  eX  =  0, 

where  the  matrices  A  and  X  are  given  in  (3. 5),  (3. 6). 

Proof :  Multiplying  (5.30)  with  IJ^Pl  and  I^Pr  leads  to 

HIPlU  +  llPRV)t  =  -{IIQlAU  +  IIQrAV)  +  tUlRiXU  +  IIRrXV) 

+(^i  -  ^r){Un  -  Vo)  +  {al  -  al){VUN  -  Wo) 

(5.33)  +2(C/,  F)p,  +  2(y,  F)p^  + 

where  BT  includes  the  boundary  terms  at  x  =  ±1.  To  obtain  (5.33)  we  have  made  use  of  Lemma  1. 

The  inviscid  terms  can  be  written 

(5.34)  IIQlAU  +  IIQrA  =  - {AU)o  +  {AU)n  -  {AV)o  +  iAV),n ■ 

Next,  we  consider  the  viscous  terms.  Both  R  =  QP  and  R  =  -f  D)5  lead  to 

IIQlPl  "QlXU  +  IIQrPp^QrXV  =  -  {X'DU)o  +  (XVU)n 

(5.35)  -{XW)o  +  {XW)m. 

By  inserting  (5.34)  and  (5.35)  into  (5.33),  neglecting  the  boundary  terms  at  x  =  ±1,  letting  F  =  0,  and  ap¬ 
plying  condition  (5.32),  we  obtain  {I^PlU -\-l^PRV)t  =  0;  i.e.,  the  approximation  (5.30)  is  conservative.  □ 

5.2.2.  Stability.  We  start  with  the  following  observation. 

Remark.  Stability  of  the  one  domain  problem  does  not  imply  stability  of  the  multiple  domain  problem. 
Stability  means  that  the  solution  can  be  estimated  in  terms  of  the  (bounded)  boundary  data.  In  a  multiple 
domain  problem,  the  boundary  data  are  made  up  of  the  solution(s)  in  the  other  domain(s).  Boundedness  of 
the  data  would  require  an  a  priori  assumption. 

The  main  result  of  this  paper  is  given  below. 

Theorem  6.  The  approximation  (5.30)  of  the  problem  (4-1)  is  both  strictly  and  strongly  stable  if 

(5.36)  al  =  cxeX,  ai  =  \{A- ptX  -  5h),  +  il±^,  5>Q, 

2  2aR  2aL 
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and  if  Theorem  3  and  5  hold.  aL^o^R  denotes  the  minimal  eigenvalue  of  P  if  R  =  QP  and  the  minimal 
eigenvalue  of  {M  +  M^)/2  if  R  =  {-S^M  -f  D)S.  The  matrices  A,  X  are  given  in  (3. 5), (3. 6). 

Proof  :  Strict  and  strong  stability  of  (5.30)  follows  if  the  interface  treatment  at  a:  =  0  is  of  a  dissipative 
nature.  For  that  reason  we  neglect  the  terms  at  the  boundaries  x  —  ±1  and  use  F  —  0.  The  energy  method 
leads  to 

JtiWUfp,  +  WVfpJ  =  -U'^{IQl  +  Qll)U  -  v^{KQr  +  Qll)V 

+eU'^{XRL  +  RlX)U  +  eV'^iXRR  +  RlX)V 
+2Ul{ai{Un  -  Vb)  +  {VUn  -  Wo)) 

(5.37)  +2Vj'{al,{Vo  -  Un)  +  (Wo  -  VUn)). 


Equations  (5. 11), (5. 12)  and  (5.13)  lead  to 

(5.38)  -U^{AQl  +  QlA)U  =  AUo  -  U^AU„ 

(5.39)  -V^{AQr  +  QlA)V  =  AVo  -  V^AVm 

(5.40)  U‘^{XRl  +  RlX)U  <  -2UoXVUo  +  2UnXVUn  -  2aLVU^XVUn 

(5.41)  V^{XRr  +  RlX)V  <  -2VoXVVo  +  2VmXVVm  -  2aRVV^ XVVo- 


By  inserting  (5.38)-(5.41)  into  (5.37)  and  neglecting  boundary  terms  at  x  =  ±1  we  obtain 

j^{\mk  +  \\y\\k)<w^^w, 


where 


W  = 


(  Un  \ 
Vo 
VU, 

V  VVo  J 


,E  = 


(  2cr[  —  A 
al  +  eX 

V  - 


-{a[  +  c^r) 
2c7^  +  A 


a 


^  -tX 


R 


^Y  +  eX 

-<^R 

0 


—a 

.■V 


\ 


aY-eX 
0 

—2aReX  J 


The  problem  (5.30)  is  strictly  and  strongly  stable  if  E  is  negative  semidefinite.  E  is  an  almost  full 
matrix;  to  obtain  explicit  stability  conditions,  simplifications  of  E  are  necessary.  The  energy  method  applied 


to  the  continuous  multiple  domain  problem  leads  to  (4.11)  which  suggests  that  the  variables 


Un- Vo  \ 

f  +I 

-I 

0 

0  ^ 

1 

Un  +  Vo 

o  1 

+I 

+I 

0 

0 

w  =  sw=-^ 
V2 

VUn  -  VVo 

0 

0 

+/ 

-I 

VUn  +  VVo  ) 

0 

0 

+/ 

are  of  interest.  The  use  of  these  variables  and  the  conservation  conditions  in  Theorem  5  leads  to 

2al  +  eX  eX  \ 

0  0 
-{aL  +  aR)eX  {aR-aL)e.X 

(aR  —  aL)€X  —(aL  +  o:R)eX  J 


El  =  SES'^  = 


f  2{2oi-A)  0 
0  0 
2(rl  +  eX  0 
\  eX  Q 


To  show  that  Ei  is  negative  semidefinite,  first  assume  that  the  first  condition  in  (5.36)  holds.  Secondly, 
add  and  subtract  the  matrix  -2(3eX  to  the  upper  left  block  in  Ei.  The  condition  for  negative  semidefiniteness 
(see  Lemma  1)  becomes 

(5.42)  yT{2{2ai  -  A)  +  2peX)yi  +  €[{Y  0  R)'^y2f[AE2  ®  B][(^  ®  Rfv^]  <  0, 
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where  B  =  R^XR,  Ae,  =  Y^EiY  and 


/  -2/? 


E2  = 


^  (l  +  2cr) 


(1  +  2(7) 
-{at  +  olr) 

OLR  -  OtL 


OLR  -  OCR 
-{oll  +  aR) 


The  first  term  in  (5.42)  is  nonpositive  if  the  seconiJ  relation  in  (5.36)  holds.  Negative  definiteness  that  implies 
<  0  is  obtained  if  the  third  relation  in  (5.36)  holds.  □ 


5.2.3.  Accuracy.  In  this  section  we  will  consider  the  accuracy  close  to  the  interface.  The  procedure  is 
similar  to  the  one  used  in  section  5.1.2  for  the  single  domain  problem.  The  problem  describing  the  deviations 
Uj  =  U{xj^t)  —  Uj{t)  and  Vj  =  V{xj,t)  —  Vj{t)  between  the  exact  continuous  solutions  and  the  discrete 
approximations  given  by  (5.30)  is 


(5.43) 


Ut  +  AVlU 

= 

eXVlU  +  Tl 

Pl  -  Vo)  +  al{{VLU)n  -  {VRV)o))eL 

C/(0) 

0 

Vt  +  AVnY 

= 

€XVIVYTr 

+ 

Pr\<^W  -U)  +  al{{VRV)o  -  {VLU)n))eR 

V{0) 

0. 

For  simplicity,  we  have  used  the  notation  U  =  U  and  V  =  V.  Note  also  that  the  terms  at  the  boundaries 
X  ~  ±1  are  neglected.  The  treatment  at  the  boundaries  x  =  ±1  has  been  discussed  in  section  5.1.2. 

Tl  and  Tr  are  the  truncation  errors  from  the  approximation  of  the  differential  operator  and  the  interface 
conditions.  The  truncation  errors  have  the  general  structure 


(  '■  ] 

f  0(Ax^"  \ 

Tl  = 

0{Ax^) 

,  Tr  = 

0{Axl) 

[  o{Ax<ir  y 

\  '■  / 

The  discussion  in  section  5.1.2  on  the  size  of  the  truncation  error  is  applicable  also  for  the  interface  problem. 
Following  the  procedure  in  section  5.1.2,  see  [31], [32],  one  splits  up  the  errors  in  two  parts,  the  first 
part  (T2,T^)  contains  the  truncation  error  of  the  internal  scheme,  and  the  second  part  contains  a  boundary 
contribution  (T|,T|)  with  one  order  lower  accuracy.  The  structure  of  these  errors  are 


(  , 

1  1 

/  : 

\ 

Tl  = 

O(Ax^) 

II 

0 

1  0 

) 

f  ' 

(  0{Ax^r 

(  «  ^ 

/ 

0{Axl) 

,  n= 

[  :  ^ 

\ 

f  0{Ax^^  '))  \ 
0 

/ 


Also  the  error  is  divided  into  two  parts;  i.e,  we  consider  U  =  -\-U^  and  V  = 

By  using  the  energy  method,  and  will  be  bounded  by  and  T^.  This  procedure  is  straightfor¬ 
ward,  entirely  similar  to  the  one  in  section  5.1.2  and  will  therefore  not  be  repeated  here.  Suffice  it  to  say 
that  the  stability  conditions  given  in  Theorem  6  lead  to 


\\U^P,  +  III^^IIph  <  OiAxl)  +  0{AxS). 
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To  bound  U'^  and  in  terms  of  T|  and  requires  use  of  the  Laplace-transform  technique.  That  analysis 
is  given  in  detail  below. 

Also  in  this  case,  we  keep  the  algebraic  complexity  down  by  considering  the  inviscid  (e  -  0)  Euler 
equations  approximated  with  the  second  order  accurate  approximation  given  by  (A.l)  and  (A.2)  in  appendix 
A.  The  problem  for  U'^  =  U  and  =  V  obtained  by  Laplace-transforming  (5.43)  becomes 


(5.44) 


SlUu  +  A([jn  —  Un  l) 

sflVo+A(yi-yo) 

slUj  +  A{Uj+i-Uj  i)/2 
SRVJ  +  K{VJ+^~Vj  x)/2 

Uj 

Vi 


2cr[(U„  -  Vo)  +  Ax^gi 
2crji(Vo  -  Un)  +  AxRgR 
0,  j  <  n  -  1 
0.  j>l, 

0,  j  ->  -00 

0,  i  00 


where  sr  =  sAxr,  sr  =  sAxr,  gL  =  0{Ax^^  ^^)  and  gR  =  0(Ax^  ^^). 
The  last  four  equations  in  (5.44)  lead  to 


(1\ 

(^\  1 

(5.45) 

II 

0 

(4)^' 

n  1  -2 

+  <^L 

1  (4)^  ”+4 

0 

\0  ) 

[0/  ' 

^  1  J 

.3\j  n 


(4) 


(5.46) 

II 

/  1  ' 

0 

1  (i^r)^  +  4 

1  (4)^  +  4 

/  0  N 
0 

1 

lo, 

/ 

[oj 

1  1  / 

-(s/4)  -  41  +  (s/Ai)2  ,  Aj  >  0 

(5.47)  4=0  ,Xj=0  , 

-is/Xj)  +  41  +  (s/4)"  .4<o 
“(s/4)  +  V  (s/4)^  ,Aj>o 

(5.48)  4=0  ,  Aj  =  0  , 

-(s/4)  -  41  +  (5/Aj)2  ,Xj  <0 

where  the  branch  of  the  square  root  is  the  one  with  a  positive  real  part  for  i?e(s)  >  0.  Also,  in  the  case  in 
which  Xj  =  0  presents  no  problem,  only  the  number  of  equations  in  (5.6)  is  reduced.  We  assume  Xj  ^  0  in 
the  following. 

The  coefficients  <tl  =  (4. 4>  4)^  s,nd  ctr  =  (aj^,  4. 4)^  Oe  determined  by  the  first  two  equations 
in  (5.44).  They,  together  with  the  first  condition  in  Theorem  5,  lead  to 


(5.49) 


E{sl,sr) 


SRh  -  Akj^  ^  +  a  -  2<7|,  24 

2(4  —  A)  srIs  +  Akr  +  A 
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A  nonsingular  E{sl^sb)  leads  via  the  Kreiss  condition  and  Parseval’s  relation;  (see  section  5.1.2)  to  the 
estimate 


WU^Wp,  <  0{Axl)  +  0{Ax%),  <  o(Axl)  +  0{Axl). 


It  still  must  be  shown  that  (5.29)  for  E,  defined  in  (5.49),  holds.  A  direct  calculation  using  (5.49), (5.47), 
and  (5.48)  leads  to 

Det{E)  =  llGj,  Gj  =  |A,f  (1  +  +  + 

3=1 

+  (sL/Aj)^y^l  +  (5ie/Aj)2. 

Let  +  (sl/Aj)^  =  r/L  +  and  +  i^R  where  7]L,r}R  are  non-negative.  A  simple 

algebraic  test  reveals  that  the  imaginary  part  and  the  real  part  of  Gj  cannot  be  zero  at  the  same  time  if 
the  inviscid  condition  for  stability  A  —  2(j£  <  0  in  Theorem  6  holds.  We  can  summarize  the  result  in  the 
following  Theorem. 

Theorem  7.  The  approximation  (5.30)  of  the  problem  (4’1)  with  e  =  0  is  second  order  accurate;  i.e., 

l|C/||p.  +  ||I^llp.<0(A:ri)  +  (9(A4), 


if  Theorem  6  holds  and  the  first  derivative  operator  V  =  P  is  given  by  (A.l)  and  (A. 2)  in  appendix  A. 

Remark.  Also  in  the  interface  case,  see  section  5.1.2,  the  procedure  to  prove  accuracy,  which  was 
exemplified  above  in  the  second  order  accurate  case,  is  general.  The  last  step  in  which  one  uses  the  Laplace- 
transform  technique  to  estimate  the  errors  U'^andV'^  is  not  necessary  i)  if  the  stencils  adjacent  to  the  interface 
have  the  same  order  of  accuracy  as  the  internal  stencil,  and  ii)  if  the  interface  conditions  are  one  order  more 
accurate. 


5.2.4.  The  discrete  multiple  domain  problem  in  conservation  form.  The  discrete  multiple 
domain  problem  (5.30)  can  be  transformed  to  conservative  form  by  multiplying  the  equations  with  In^i  0 
T{RS)  ^,/m+i  ^T{RS)  respectively.  The  result  is 


(5.50) 


Ut  +  Pr,\QLF^  -eRLF'^)  = 
Vt  +  Pp\QRF^  -eRRF^)  = 


+  il/2)Pj^^[{Fl  -  F^) 

+  (l  +  2c7)e(Fif-Fj)-Ff] 

-  (l/2)P^i[(Fj-Ff) 

-  il  +  2cr)eiFX-Fl)  +  F^], 


where  F'^  =  F^  -  eF^  and 


Pf  =  (Sh  +  epfBT  ^){Un  -  Vo),  Fi  =  (Sis  +  epfSf  i)(yo  -  t4). 

In  (5.50),  the  forcing  terms  and  the  boundary  conditions  at  a;  ±  1  are  neglected. 

6.  Numerical  experiments.  By  making  one-dimensional  computations  using  the  nonlinear  Euler  and 
Navier- Stokes  equations,  we  can  check  whether  the  theoretical  conclusions  drawn  from  the  analysis  of  the 
constant  coefficient  problem  agree  with  the  results  obtained  in  practice. 

In  the  calculations  below,  we  use  the  second  order  scheme  (given  in  (A.1)-(A.5))  and  the  fourth  and 
sixth  order  schemes  reported  in  [24].  To  integrate  in  time,  a  five-stage  fourth-order  RK  scheme  [33]  has  been 
used.  Consider  the  stability  condition  (5.36).  In  the  calculations  below  we  have  used  cr  =  —1/2  and  the 
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Fig.  6.1.  L2  Errors  in  calculations  using  the  Euler  equations. 

conservative  estimate  cr£  =  A/2  -  51  where  5  is  determined  through  tests.  Often  we  use  5  —  1.0.  Equation 
(5.32)  has  been  used  to  determine  the  other  parameters. 

First,  we  consider  a  sound  propagation  problem.  The  computational  results,  obtained  using  the  nonlinear 
Euler  equations  at  Mach  number  0.5,  are  compared  with  an  exact  solution  of  the  linearized  problem.  In 
Figure  6.1  The  errors  for  second,  fourth,  and  sixth  order  schemes  using  one  domain  (IDom),  four  uniform 
domains  (4Dom),  and  eight  randomly  spaced  domains  (Rand)  are  shown.  Clearly,  the  order  of  accuracy  is 
independent  of  the  presence  and  location  of  the  interfaces.  Due  to  the  small  amplitudes  (oc  10  ’’)  used  in 
the  sixth  order  cases,  we  encounter  round  off,  which  can  be  seen  as  the  kink  on  the  sixth  order  results. 

Next,  we  consider  a  viscous  shock  propagation  problem  at  Mach  number  2.0  and  Reynolds  number  150. 
The  exact  solution  of  the  Navier-Stokes  equation  for  this  case  can  be  found  in  [34].  In  Figure  6.2,  the  errors 
for  second,  fourth,  and  sixth  order  schemes  using  eight  uniform  domains  (Unif)  and  eight  randomly  spaced 
domains  (NonU)  are  shown.  Also  in  this  case,  the  order  of  accuracy  is  independent  of  the  location  of  the 
interfaces. 

The  curves  in  the  sixth  order  case  are  not  straight,  see  Figure  6.2.  The  curves  are  formed  as  a  mean 
value  of  15  simulations  where  different  wave  speeds  ws  from  —0.25  to  0.5  are  used.  The  individual  results 
for  each  wave  speed  are  given  in  Tables  6.1  and  6.2.  Note:  ws  —  0  is  stationary  shock,  and  we  have  subsonic 
wave  speeds  fox  ws  <  0.3.  The  results  from  uniform  grid  calculations  are  shown  in  Table  6.1;  results  from 
nonuniform  grid  calculations  are  given  in  Table  6.2.  The  convergence  rate  between  the  two  grids  is  listed. 
The  asymptotic  limit  approaches  -6.  Note  that  the  trends  are  identical  between  the  nonuniform  and  uniform 
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Viscous  Shock  Propagation 


Fig.  6.2.  L2  Errors  in  calculations  using  the  Navier- Stokes  equations. 


cases. 

In  Figure  6.3,  the  propagating  shock  (ws  =  0.25)  for  four  different  times  is  shown.  In  this  case,  the  sixth 
order  scheme  and  24  gridpoints  were  used  in  eaeh  domain. 

Finally  we  will  discuss  two  additional  questions  concerning  accuracy  and  stability/efBciency.  To  investi¬ 
gate  the  influence  of  interface  conditions  on  accuracy,  we  made  the  calculations  illustrated  in  Table  6.3.  The 
calculations  are  run  to  a  physical  time  T  =  3  at  Mach  number  2.0  and  Reynolds  number  Re  =  250.  The 
sixth  order  SBP  scheme  is  used,  and  the  number  of  total  points  is  289  evenly  distributed  on  the  interval 
-1/2  <  X  <  1.  The  parameter  in  the  study  is  the  number  of  subdomains,  keeping  the  total  number  of 
intervals  constant.  The  number  of  subdomains  ranges  from  1  to  24.  For  the  case  of  24  subdomains,  the 
spatial  operator  involves  12  boundary  stencils  (fifth-order)  and  one  sixth-order  interior  stencil.  No  further 
divisions  are  possible  when  using  the  sixth-order  SBP  operator.  Note  that  this  case  is  only  marginally  less 
accurate  than  the  single  domain  case,  for  which  the  most  points  are  discretized  with  sixth-order  stencils. 

The  previous  study  indicates  that  there  is  little  loss  of  accuracy  when  subdividing  the  domain.  There 
are,  however,  other  costs  associated  with  domain  subdivision.  Introduction  of  additional  interfaces  into  the 
domain  changes  the  resulting  eigenspectrum  of  the  semidiscrete  operator.  In  [22],  a  reduction  in  the  effective 
CFL,  when  using  a  penalty  boundary  procedure,  was  observed.  We  experience  a  similar  reduction  in  the 
stability  envelop  as  the  number  of  subdomains  is  increased. 

In  Table  6.4,  a  study  compares  the  effective  CFL  of  a  singledomain  calculation,  with  those  from  a 
comparable  grid  divided  into  eight  subdomains.  Plotted  are  the  errors,  and  the  maximum  stable  CFL  as  a 
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wave 

speed 

-0.2500 

-0.2000 

-0.1500 

-0.1000 

-0.0500 

0.0000 

0.0500 

0.1000 

0.1500 

0.2000 

0.2500 

0.3000 

0.3500 

0.4000 

0.4500 

0.5000 


a:  96 

b:128 

-4.4460 

-3.1743 

-1.5344 

-3.4447 

-4.7040 

-3.3093 

-5.6203 

-2.0256 

-4.9470 

-7.5646 

-4.7062 

-5.9644 

-4.9922 

-3.8538 

-0.7633 

-4.8735 


-4.4858  -4.8856 
-4.3626  -4.6413 

-6.5479  -2.8759 

-5.0035  -6.5253 

-4.8410  -5.0488 

-6.0257  -5.0116 
-3.0687  -3.8456 

-6.6065  -6.2051 

-5.1708  -5.5378 
-5.5734  -6.0670 

-3.0715  -3.4963 

-6.8890  -6.1815 

-5.0159  -5.5773 

-5.9798  -5.3015 

-2.4389  -5.8286 

-5.9665  -5.7257 


256 

384 

-5.2375 

-4.7217 

-7.4052 

-5.6487 

-5.4075 

-5.4064 

-6.8111 

-5.6289 

-5.4018 

-5.6418 

-4.1048 

-6.1886 

-5.4589 

-6.3515 

-4.9988 

-5.8033 


384 

512 

-5.5610 

-5.2825 

-5.6489 

-6.2999 

-5.6641 

-5.6608 

-6.4191 

-5.2130 

-5.7503 

-5.8898 

-5.8041 

-6.2484 

-5.1898 

-5.8163 

-7.4170 

-5.8160 


Table  6.1 

UNIFORM  grid,  8  subdomains,  refinements  between  grids  a:b,  6th  order  explicit;  CFL  =  O.S. 


NONUNIFORM  g 
max/min  ratio  is  6.47. 


wave  a:  96  128  192  2DD  884 

speed  b:128  192  256  384  |  512 

-0.2500  I  -3.9508  I  -5.0473  I  -5.6119  I  -5.9785  |  -4.5714 

-0.2000  -3.8816  -4.6038  -4.8316  -6.3334  -6.9001 

-0.1500  -4.4373  -4.8932  -5.1011  -5.4514  -5.8566 

-0.1000  -7.6431  -3.4782  -7.2908  -5.8925  -5.9310 

-0.0500  -3.2678  -7.4661  -3.7127  -6.8562  -6.2791 

0.0000  -4.8406  -4.7810  -5.2543  -5.5640  -5.7606 

0.0500  -7.8667  -2.9978  -7.8219  -4.5670  -4.7133 

0.1000  -4.2532  0.6385  -2.8840  -3.9126  -6.1144 

0.1500  -1.4577  -5.1589  -6.7711  -5.0271  -5.2246 

0.2000  -4.2245  -4.6373  -4.5406  -5.2829  -5.9054 

0.2500  -2.9734  -4.0155  -5.4352  -5.4569  -5.6145 

0.3000  -4.2383  -3.2435  -4.3861  -5.5030  -5.8943 

0.3500  -4.1902  -3.5835  -3.6667  -4.6596  -5.5589 

0.4000  -3.4505  -3.4125  -4.8703  -5.3783  -6.3374 

0.4500  -2.7380  -2.9597  -3.6868  -4.5673  -4.8980 

0.5000  -4.1279  -3.7378  -3.8465  -4.8352  -6.1756 

Table  6.2 

8  subdomains  generated  randomly f  refinements  between  grids  a:b,  6th  order  explicit,  CFL 


viscous  Shock  Propagation 


Fig.  6.3.  Viscous  shock  propagation,  a  domain  with  randomly  spaced  interfaces. 


Subdomains 

LOGio^rror 

1 

-4.527 

2 

-4.584 

4 

-4,457 

8 

-4.643 

12 

-4.313 

16 

-4.467 

18 

-4.342 

24 

-4.358 

Table  6.3 


Variation  of  L2  error  on  number  of  subdomains  with  grid  density  constant. 


function  of  Reynolds  number  for  the  two  cases.  Note  that  while  the  errors  are  nearly  equivalent  for  the  two 
test  cases,  the  maximum  CFL  for  the  single  domain  case  is  nearly  a  factor  of  two  larger. 

7.  Summary  and  conclusions.  We  have  analyzed  boundary  conditions  and  interface  conditions  for 
the  one-dimensional  Euler  and  Navier-Stokes  equations.  Both  the  continuous  and  semi-discrete  problems 
have  been  considered. 

We  have  considered  summation-by-parts  operators  and  derived  strictly  and  strongly  stable  boundary 
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Re 

LOGi^error 

OF  L-fnax 

LOGioerror 

GF  Lqjiax 

1000 

-2.154 

0.55 

DNC 

900 

-2.242 

0.55 

-2.265 

0.30 

800 

-2.347 

0.55 

-2.376 

0.30 

700 

-2.477 

0.60 

-2.517 

0.30 

600 

-2.637 

0.60 

-2.698 

0.30 

500 

-2.841 

0.60 

-2.935 

0.30 

300 

-3.429 

0.65 

-3.617 

0.30 

200 

-4.027 

0.65 

-4.185 

0.35 

100 

-5.741 

0.60 

-5.699 

0.35 

40 

-7.892 

0.50 

-7.331 

0.20 

20 

-9.535 

0.45 

-8.637 

0.20 

10 

-10.968 

0.40 

-10.665 

0.18 

Table  6.4 

Variation  of  CFL  number  and  L2  error  with  Reynolds  number  for  single  and  multiple  domain  cases. 


and  interface  conditions  for  the  Euler  and  Navier-Stokes  equations.  We  have  also  considered  the  question 
of  accuracy,  both  in  the  general  case  and  more  specifically  for  a  second  order  accurate  approximation  of  the 
Euler  equations. 

The  interface  conditions  are  stable  and  conservative  even  if  the  finite  difference  operators  and  mesh  sizes 
vary  from  domain  to  domain.  Numerical  experiments  which  include  a  sound  propagating  problem  and  a 
viscous  shock  propagating  problem  show  that  the  new  conditions  lead  to  accurate  and  stable  results  for  the 
corresponding  nonlinear  problems  also. 

It  was  also  shown  by  numerical  experiments  that  there  is  little  loss  of  accuracy  associated  with  domain 
subdivision.  However,  the  introduction  of  interfaces  into  the  domain  changed  the  eigenspectrum  of  the 
semidiscrete  operator  and  caused  a  reduction  of  the  CFL  number  by  approximately  a  factor  of  two. 


Appendix  A.  Stencils.  We  now  present  a  few  examples  of  the  specific  form  of  the  stencils  that  have 
the  SBP  property.  For  more  accurate  stencils,  see  [24].  The  second  order  accurate  discretization  matrix  that 
approximates  the  first  derivative  P  =  P  is 


(A.1) 


-2  2 

-1  0  1 


-1  0  1 
-2  2 
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where 


(A.2) 


1 - 

1 _ 

■-1  1 

1 

-1  0  1 

.0=1 

.  . 

1 

-1  0  1 

1 

2  - 

- 1 

rH 

1 

_ 1 

The  second  order  accurate  discretization  matrix  that  approximates  the  second  derivative  —  F 


D)S  is 


(A.3) 


1  -2  1 

1  --2  1 


1-2  1 
1  -2  1 


where 


(A.4)  5 


and 


1  «_2  2 

9  9  9 

_2  10  _10 

9  9  9 

2  _10  19  _i 

9  9  9-^ 

-1  2  -1 


-1  2  -1 

1  19  _2 

•^9  9  9 

_  10  _  10  2 

9  9  9 

_  2  _  2  4 

9  9  9 


(A,5) 


Aa: 


The  matrix  M  can  be  shown  to  be  positive  definite  (and  symmetric). 
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